Prep
setwd('/project/CRI/DeBerardinis_lab/lcai/NSCLC_Metabolism/scripts/final/')
.libPaths(c(.libPaths(),'~/R/x86_64-unknown-linux-gnu-library/3.2/'))
# libraries
library(GGally)
library(openxlsx)
library(Hmisc)
library(corrplot)
library(RColorBrewer)
library(GSVA)
library(ppcor)
library(reshape2)
library(gtable)
library(grid)
library(gridExtra)
library(factoextra)
library(ggrepel)
library(heatmap3)
library(stringr)
library(plotrix)
library(mclust)
library(MASS)
library(sjPlot)
library(ComplexHeatmap)
library(ClassComparison)
library(pracma)
library(ggpubr)
library(cowplot)
# load data
c2cgp<-readRDS('data/c2cgp.rds')
c2cp<-readRDS('data/c2cp.rds')
dat<-readRDS('data/metabolic_profiling_data.rds')
ge<-readRDS('data/illumina_gene_level.rds')
common<-intersect(rownames(ge),rownames(dat))
ge<-ge[common,]
# ssGSEA.results<-gsva(expr = t(ge),gset.idx.list = c2cgp,method="ssgsea",rnaseq=F)
# The line above is commented because it takes too long to run, the pre-calculated results were saved and loaded below
ssGSEA.results<-readRDS('data/ssGSEA_result.rds')
sig.ne <- read.xlsx('data/Zhang_2018_NE_signature.xlsx')
pc<-readRDS('data/34glucose_PC.rds')
mutations<-readRDS('data/mutations.rds')
dat.invivo<-readRDS('data/invivo_data.rds')
mut<-readRDS('data/invivo_mutation.rds')
lkb1<-readRDS('data/LKB1_OE.rds')
S4.data<-readRDS('data/Table_S4.rds')
AUC<-read.xlsx('data/table6_chemical_screen_data.xlsx',sheet = 'AUC')
emt.signature<-read.csv('data/EMT_signature.csv')
emt.signature<-unique(emt.signature$Gene.Symbol)
emt.signature<-gsub(' ','',fixed = T,emt.signature)
pem<-readRDS('data/pem_validation.rds')
wb.misty<-loadWorkbook('data/Misty_xenograft.xlsx')
wb<-loadWorkbook('data/All Metabolic data - 2.16.2015R1.xlsx')
CDH1.data<-readRDS('data/CDH1_multiomics.rds')
load('data/additional_serine_data.RData')
col2<-colorRampPalette(rev(brewer.pal(9,'RdBu')))(100)
# Shared functions
cor_fun <- function(data, mapping, method="pearson", ndp=2, sz=5, stars=TRUE, ...){
xVal<-mapping_string(mapping$x)
yVal<-mapping_string(mapping$y)
data <- na.omit(data[,c(xVal,yVal)])
x<-data[,xVal]
y<-data[,yVal]
corr <- cor.test(x, y, method=method)
est <- corr$estimate
if(stars){
stars <- c("***", "**", "*", "")[findInterval(corr$p.value, c(0, 0.001, 0.01, 0.05, 1))]
lbl <- paste0(round(est, ndp), stars)
}else{
lbl <- round(est, ndp)
}
ggplot(data=data, mapping=mapping) +
annotate("text", x=mean(x), y=mean(y), label=lbl)+
theme(panel.grid = element_blank())
}
col.grad<-function(x){
scales::alpha(colorRampPalette(rev(brewer.pal(5,'RdBu')))(length(x)),alpha = 0.5)[rank(x)]
}
colorlegend <- function(
colbar,
labels,
at = NULL,
xlim = c(0, 1),
ylim = c(0, 1),
vertical = TRUE,
ratio.colbar = 0.4,
lim.segment = "auto", # NOTE: NULL treated as "auto"
align = c("c", "l", "r"),
addlabels = TRUE,
...)
{
if (is.null(at) && addlabels) {
at <- seq(0L, 1L, length = length(labels))
}
if (is.null(lim.segment) || lim.segment == "auto") {
lim.segment <- ratio.colbar + c(0, ratio.colbar * .2)
}
if (any(at < 0L) || any(at > 1L)) {
stop("at should be between 0 and 1")
}
if (length(lim.segment) != 2) {
stop("lim.segment should be a vector of length 2")
}
if (any(lim.segment < 0L) || any(lim.segment > 1L)) {
stop("lim.segment should be between 0 and 1")
}
align <- match.arg(align)
xgap <- diff(xlim)
ygap <- diff(ylim)
len <- length(colbar)
rat1 <- ratio.colbar
rat2 <- lim.segment
if (vertical) {
at <- at * ygap + ylim[1]
yyy <- seq(ylim[1], ylim[2], length = len + 1)
rect(rep(xlim[1], len), yyy[1:len],
rep(xlim[1] + xgap * rat1, len), yyy[-1],
col = colbar, border = colbar)
rect(xlim[1], ylim[1], xlim[1] + xgap * rat1, ylim[2], border = "black")
segments(xlim[1] + xgap * rat2[1], at, xlim[1] + xgap * rat2[2], at)
if (addlabels) {
pos.xlabel <- rep(xlim[1] + xgap * max(rat2, rat1), length(at))
switch(align,
l = text(pos.xlabel, y = at, labels = labels, pos = 4, ...),
r = text(xlim[2], y = at, labels = labels, pos = 2, ...),
c = text((pos.xlabel + xlim[2]) / 2, y = at, labels = labels, ...),
stop("programming error - should not have reached this line!")
)
}
} else {
at <- at * xgap + xlim[1]
xxx <- seq(xlim[1], xlim[2], length = len + 1)
rect(xxx[1:len], rep(ylim[2] - rat1 * ygap, len),
xxx[-1], rep(ylim[2], len),
col = colbar, border = colbar)
rect(xlim[1], ylim[2] - rat1 * ygap, xlim[2], ylim[2], border = "black")
segments(at, ylim[2] - ygap * rat2[1], at, ylim[2] - ygap * rat2[2])
if (addlabels) {
pos.ylabel <- rep(ylim[2] - ygap * max(rat2, rat1), length(at))
switch(align,
l = text(x = at, y = pos.ylabel, labels = labels, pos = 1, ...),
r = text(x = at, y = ylim[1], labels = labels, pos = 2, ...),
c = text(x = at, y = (pos.ylabel + ylim[1]) / 2, labels = labels, ...),
stop("programming error - should not have reached this line!")
)
}
}
}
addLegend<-function(x,name){
colbar<-col.grad(sort(x))
labels<-cl.lim<-round(range(x,na.rm = T),1);cl.length<-100
cl.align.text = "c"
par(mar=c(0.5,3.2,2,3.2),xpd=T)
plot.new()
title(name,line = 0)
colorlegend(colbar = colbar, labels = labels,
offset = 1, ratio.colbar = 0.5, cex = 1,
vertical = F,
align = cl.align.text)
}
tri_plot<-function(x){
addLegend(x[,3],name=names(x)[3])
par(mar=c(3,3,0.2,3))
plot(x[,1],x[,2],col=col.grad(x[,3]),pch=19,xlab=names(x)[1],ylab=names(x)[2],mgp=c(2,1,0),oma=c(2,0,2,0))
addLegend(x[,1],name=names(x)[1])
par(mar=c(3,3,0.2,3))
plot(x[,2],x[,3],col=col.grad(x[,1]),pch=19,xlab=names(x)[2],ylab=names(x)[3],mgp=c(2,1,0),oma=c(2,0,2,0))
addLegend(x[,2],name=names(x)[2])
par(mar=c(3,3,0.2,3))
plot(x[,3],x[,1],col=col.grad(x[,2]),pch=19,xlab=names(x)[3],ylab=names(x)[1],mgp=c(2,1,0),oma=c(2,0,2,0))
}
Isotopologue.Distribution<-function(Key,Name){
dat.key<-dat[,grep(Key,grep("m[0-9]$",names(dat),value = T))]
dat.key<-setNames(melt(dat.key[complete.cases(dat.key),]),c('Isotopologue','% Enrichment'))
dat.key$tracer<-factor(ifelse(grepl('Q',dat.key$Isotopologue),'Glutamine','Glucose'),levels=c('Glucose','Glutamine'))
dat.key$timepoint<-factor(ifelse(grepl('24',dat.key$Isotopologue),'24h','6h'),levels=c('6h','24h'))
dat.key$Isotopologue<-sapply(dat.key$Isotopologue,FUN = function(x){
tmp<-unlist(strsplit(split = 'm',as.character(x)))
return(paste0('m+',tmp[length(tmp)]))
})
ggplot(dat.key, aes(x=Isotopologue, y=`% Enrichment`)) +
geom_violin(trim = T,fill='#6e016b',alpha=0.5)+
geom_boxplot(width=0.1,outlier.shape = NA)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
theme_bw()+
facet_grid(tracer~timepoint)+
ggtitle(paste0(Name,' Labeling'))
}
Fig 1
1B
ggpairs(dat[,c('Glc','Lac','Gln','Glu')],
upper=list(continuous=cor_fun),
lower=list(continuous = wrap("points", alpha = 0.5,colour='#2c7fb8')))+
theme_bw()+theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
1C-F
LacGlc_ssGSEA<-merge(dat[,'Lac/Glc',drop=F],t(ssGSEA.results),by = 0,all = F)
LacGlc_ssGSEA.r<-apply(LacGlc_ssGSEA[,-c(1:2)],2,FUN=function(x) cor.test(x,LacGlc_ssGSEA$`Lac/Glc`)$estimate)
lacGlc_hypoxia<-LacGlc_ssGSEA.r[grep('HYPOXIA',names(LacGlc_ssGSEA.r))]
lacGlc_hypoxia<-lacGlc_hypoxia[grep('DN',names(lacGlc_hypoxia),invert = T)]
layout(matrix(c(1,2,4,3,2,4),byrow = T,ncol = 3,nrow = 2),widths=c(3,2,2), heights=c(1,1))
plot(density(lacGlc_hypoxia,cut = range(lacGlc_hypoxia),bw = 0.05),col='#810f7c',lwd=2,
main='Correlation between Lac/Glc and Hypoxia Gene Signature Enrichment Scores',
cex.main=0.9)
points(lacGlc_hypoxia,jitter(rep(0.55,length(lacGlc_hypoxia)),amount = 0.1),col=scales::alpha('#810f7c',0.4),pch=19)
tmp<-LacGlc_ssGSEA[,c('WINTER_HYPOXIA_UP','Lac/Glc')]
tmp.cor<-rcorr(as.matrix(tmp))
plot(tmp,
main=paste0('r = ',round(LacGlc_ssGSEA.r['WINTER_HYPOXIA_UP'],2),', pv = ',signif(tmp.cor$P[1,2],2)),
pch=19,col=scales::alpha(colour = '#d7301f',alpha = 0.5),
xlab='WINTER_HYPOXIA_UP enrichment score')
winter_hypoxia_r<-rcorr(cbind(dat[common,]$`Lac/Glc`,ge[common,intersect(c2cgp$WINTER_HYPOXIA_UP,colnames(ge))]))$r[1,-1]
plot(density(winter_hypoxia_r,bw = 0.05,cut = range(winter_hypoxia_r)),col='#810f7c',lwd=2,
main='Correlation between Lac/Glc and WINTER_HYPOXIA_UP genes',cex.main=0.9)
set3<-winter_hypoxia_r[!names(winter_hypoxia_r)%in%c(c2cp$REACTOME_GLYCOLYSIS,'LDHA')]
points(set3,jitter(rep(0.55,length(set3)),amount = 0.1),col=scales::alpha('#810f7c',0.1),pch=19)
set1<-winter_hypoxia_r[names(winter_hypoxia_r)%in%c2cp$REACTOME_GLYCOLYSIS]
points(set1,jitter(rep(0.55,length(set1)),amount = 0.1),col=scales::alpha('black',0.8),pch=19)
set2<-winter_hypoxia_r[names(winter_hypoxia_r)=='LDHA']
points(set2,jitter(rep(0.55,length(set2)),amount = 0.1),col=scales::alpha('red',0.8),pch=19)
legend('topright',legend = c('REACTOME_GLYCOLYSIS','LDHA','other'),bty='n',bg = 'transparent',cex = 0.8,pch=19,
col=scales::alpha(alpha = c(0.8,0.8,0.1),colour = c('black','red','#810f7c')))
comp.score <- function(dat, sig) {
ind <- match(sig$Symbol,rownames(dat))
ind <- ind[!is.na(ind)]
sig <- sig[sig$Symbol %in% rownames(dat), ]
score1 <- apply(dat[ind, ], 2, function(d) cor.test(sig[, 'NE.class.mean.expression'], d)$estimate)
score2 <- apply(dat[ind, ], 2, function(d) cor.test(sig[, 'Non-NE.class.mean.expression'], d)$estimate)
score <- (score1 - score2)/2
score
}
ne_score<-comp.score(t(ge),sig = sig.ne)
tmp<-data.frame(ne_score,dat[common,'Lac/Glc'])
tmp<-tmp[complete.cases(tmp),]
tmp.cor<-rcorr(as.matrix(tmp))
plot(tmp,main=paste0('r = ',round(tmp.cor$r[1,2],2),', pv = ',signif(tmp.cor$P[1,2],2)),
pch=19,col=scales::alpha(colour = '#d7301f',alpha = 0.5),
xlab='neuroendocrine score',ylab='Lac/Glc')
1G
corrplot(rcorr(as.matrix(data.frame('Glu/Gln'=dat[common,c('Glu/Gln')],ge[,c('GLS','GLS2')],check.names = F)))$r,cl.pos = 'n',
col=col2,addCoef.col = "black",diag = F,type = 'lower',method = 'color',tl.col='black')
1H
corrplot(rcorr(as.matrix(dat[,c('Glc','Lac','Lac/Glc','Gln','Glu','Glu/Gln',
'Day3/Day1','Day5/Day1',
'Day1-G','Day3-G','Day5-G','Day1-Q','Day3-Q','Day5-Q')]))$r,
col=col2,diag = F,type = 'upper',method = 'color',addCoef.col = "black",tl.col='black')
1I top
make_cormat<-function(x){
corMat<-cor(x)
pcorMat<-pcor(x)$estimate
colnames(pcorMat)<-rownames(pcorMat)<-colnames(corMat)
return(list(corMat,pcorMat))
}
plot_cor<-function(corMat,pcorMat,cex=1){
par(mfrow=c(1,2))
corrplot(corMat,col=col2,addCoef.col = "black",method = 'color',cl.pos = 'n',
type = 'upper',tl.cex = cex,tl.col = 'black',diag = F,title = 'pairwise correlation',mar=c(2,2,4,2))
corrplot(pcorMat,col=col2,addCoef.col = "black",method = 'color',cl.pos = 'n',
type = 'upper',tl.cex = cex,tl.col = 'black',diag = F,title = 'partial correlation',mar=c(2,2,4,2))
}
par(mfrow=c(1,2))
tmp<-make_cormat(na.omit(dat[,c('Lac','Glc','Day5-G')]));plot_cor(tmp[[1]],tmp[[2]])
1J top
par(mfrow=c(1,2))
tmp<-make_cormat(na.omit(dat[,c('Lac','Glc','Gln')]));plot_cor(tmp[[1]],tmp[[2]])
1I-J bottom
layout(matrix(1:12,ncol = 6,byrow = F),widths=rep(2,6), heights=c(2, 5))
tri_plot(dat[,c('Glc','Lac','Day5-G')])
tri_plot(dat[,c('Glc','Lac','Gln')])
Fig 2
2C
Isotopologue.Distribution(Key = 'Cit',Name='Citrate')
2D-E
c1<-rcorr(as.matrix(dat[,c('CitG6m0','CitQ6m5')]))
c2<-rcorr(as.matrix(dat[,c('CitG6m2','CitQ6m4')]))
marrangeGrob(list(ggplot(dat[,c('CitG6m0','CitQ6m5')],mapping = aes(x=CitG6m0,y=CitQ6m5))+
geom_point(colour='#6e016b',alpha=0.5,size=2.5)+theme_bw()+
ggtitle(paste0("Glutamine dependent reductive carboxylation\nr = ",round(c1$r[1,2],2),', pv = ',signif(c1$P[1,2],2))),
ggplot(dat[,c('CitG6m2','CitQ6m4')],mapping = aes(x=CitG6m2,y=CitQ6m4))+
geom_point(colour='#6e016b',alpha=0.5,size=2.5)+theme_bw()+
ggtitle(paste0("Glutamine dependent anaplerosis\nr = ",round(c2$r[1,2],2),', pv = ',signif(c2$P[1,2],2)))),
nrow=2, ncol=1, top=NULL)
2F
metabolite<-sapply(colnames(dat),FUN = function(x) unlist(strsplit(x,split = "m[0-9]"))[1])
select<-names(table(metabolite))[table(metabolite)>1]
dat.sub<-dat[,metabolite%in%select]
metabolite.sub<-sapply(colnames(dat.sub),FUN = function(x) unlist(strsplit(x,split = "m[0-9]"))[1])
G6<-grep('G6',unique(metabolite.sub),value = T)
Q6<-grep('Q6',unique(metabolite.sub),value = T)
G24<-grep('G24',unique(metabolite.sub),value = T)
Q24<-grep('Q24',unique(metabolite.sub),value = T)
nutrient.contribution<-lapply(list(G6,Q6,G24,Q24),FUN=function(c){
setNames(as.data.frame(do.call(cbind,lapply(c,FUN=function(x){
apply(dat.sub[,metabolite.sub==x],1,FUN=function(y){
sum((0:(length(y)-1))*(y/100))/(length(y)-1)
})
}))),nm = c)
})
names(nutrient.contribution)<-c('G6','Q6','G24','Q24')
tmp6<-setNames(data.frame(melt(nutrient.contribution$G6[,gsub('Q','G',colnames(nutrient.contribution$Q6))]),melt(nutrient.contribution$Q6)[,2]),
nm = c('experiment','[U-13C]glucose','[U-13C]glutamine'))
tmp24<-setNames(data.frame(melt(nutrient.contribution$G24[,gsub('Q','G',colnames(nutrient.contribution$Q24))]),melt(nutrient.contribution$Q24)[,2]),
nm = c('experiment','[U-13C]glucose','[U-13C]glutamine'))
tmp<-data.frame(time=rep(c('6h','24h'),c(nrow(tmp6),nrow(tmp24))),rbind(tmp6,tmp24),check.names = F)
tmp$experiment<-gsub('G',' ',tmp$experiment)
tmp$experiment<-gsub("$","h",tmp$experiment)
tmp$experiment<-factor(tmp$experiment)
tmp$experiment<-factor(tmp$experiment,levels=levels(tmp$experiment)[c(2,4,8,6,1,3,7,5)])
tmp$sum<-tmp$`[U-13C]glucose`+tmp$`[U-13C]glutamine`
ggplot(tmp,aes(x=sum))+geom_density()+facet_wrap(~experiment,nrow = 2)+theme_bw()+xlab("sum of glucose and glutamine fractional carbon contribution")
2G
h<-hclust(as.dist(1-abs(rcorr(as.matrix(dat),type = 'pearson')$r)),method = 'ward.D2')
fviz_dend(h, k = 23, k_colors = "jco",main = 'Clustering Metabolic Features',
type = "phylogenic", repel = TRUE)
Fig 4
4A
# KRAS mutations
KRAS<-mutations[,'KRAS']
KRAS.mut<-which(KRAS!="")
KRAS.wt<-which(KRAS=="")
KRAS.ms<-grep('MS',KRAS,ignore.case = F)
aa_change<-str_extract(pattern = "[pP].[ ]?[A-Z][0-9]{1,5}[A-Z]", string = KRAS[KRAS.ms])
aa_site<-as.numeric(str_extract(pattern = "[0-9]+",string = aa_change))
KRAS.ms12<-KRAS.ms[aa_site==12]
KRAS.ms13<-KRAS.ms[aa_site==13]
KRAS.ms61<-KRAS.ms[aa_site==61]
# EGFR mutations
EGFR<-mutations[,'EGFR']
EGFR.mut<-which(EGFR!="")
exon19.del<-str_extract(pattern = "[A-Z]7[0-9]+[-_][A-Z]7[0-9]+[ ]*del",string = EGFR)
EGFR.wt<-which(EGFR=='')
EGFR.exon19del<-which(!is.na(exon19.del))
# STK11 mutations
STK11<-mutations[,'STK11']
STK11.mut<-which(STK11!="")
STK11.wt<-which(STK11=="")
STK11.ms<-grep('MS',STK11,ignore.case = F)
STK11.ns<-grep('NS',STK11,ignore.case = F)
STK11.ss<-grep('SS',STK11,ignore.case = F)
STK11.fs<-grep('fs;',STK11,ignore.case = F)
STK11.ms<-setdiff(STK11.ms,STK11.fs)
KRAS_STK11<-intersect(KRAS.mut,STK11.mut)
select.mut<-do.call(cbind,lapply(list(KRAS.wt,KRAS.mut,KRAS.ms12,KRAS.ms13,KRAS.ms61,EGFR.wt,EGFR.mut,EGFR.exon19del,STK11.wt,STK11.mut,STK11.ms,STK11.ns,STK11.fs,STK11.ss), FUN=function(x) (1:nrow(mutations))%in%x))
colnames(select.mut)<-c('KRAS.wt','KRAS.mut','KRAS.ms12','KRAS.ms13','KRAS.ms61','EGFR.wt','EGFR.mut','EGFR.exon19del','STK11.wt','STK11.mut','STK11.ms','STK11.ns','STK11.fs','STK11.ss')
specific.mut.mat<-matrix(NA,ncol=ncol(dat),nrow = ncol(select.mut)-3,dimnames = list(grep('wt',colnames(select.mut),invert = T,value = T),colnames(dat)))
for(i in 1:nrow(specific.mut.mat)){
pick<-grep('wt',colnames(select.mut),invert = T,value = T)[i]
pick.gene<-unlist(strsplit(pick,split = '.',fixed = T))[1]
pick.ctrl<-paste0(pick.gene,'.wt')
for(j in 1:ncol(dat)){
specific.mut.mat[i,j]<-t.test(dat[select.mut[,pick.ctrl],j],dat[select.mut[,pick],j])$p.value
}
}
library(RColorBrewer)
KRAS<-ifelse(select.mut[,'KRAS.wt'],'white','brown')
KRAS.col<-brewer.pal(3,'Set1')
KRAS[select.mut[,'KRAS.ms12']]<-KRAS.col[1]
KRAS[select.mut[,'KRAS.ms13']]<-KRAS.col[2]
KRAS[select.mut[,'KRAS.ms61']]<-KRAS.col[3]
EGFR<-ifelse(select.mut[,'EGFR.wt'],'white','brown')
EGFR.col<-brewer.pal(3,'Set2')
EGFR[select.mut[,'EGFR.exon19del']]<-EGFR.col[1]
STK11<-ifelse(select.mut[,'STK11.wt'],'white','brown')
STK11.col<-brewer.pal(4,'Set3')
STK11[select.mut[,'STK11.ms']]<-STK11.col[1]
STK11[select.mut[,'STK11.ns']]<-STK11.col[2]
STK11[select.mut[,'STK11.fs']]<-STK11.col[3]
STK11[select.mut[,'STK11.ss']]<-STK11.col[4]
KRAS_col<-ifelse(select.mut[,'KRAS.wt'],'#f1eef6','#dd1c77')
KRAS.ms12_col<-ifelse(!select.mut[,'KRAS.ms12'],'#f1eef6',"#df65b0")
KRAS.ms13_col<-ifelse(!select.mut[,'KRAS.ms13'],'#f1eef6',"#df65b0")
KRAS.ms61_col<-ifelse(!select.mut[,'KRAS.ms61'],'#f1eef6',"#df65b0")
EGFR_col<-ifelse(select.mut[,'EGFR.wt'],'#f1eef6','#dd1c77')
EGFR.exon19del_col<-ifelse(!select.mut[,'EGFR.exon19del'],'#f1eef6','#df65b0')
STK11_col<-ifelse(select.mut[,'STK11.wt'],'#f1eef6','#dd1c77')
KRAS_STK11_col<-ifelse(!(select.mut[,'STK11.mut'] & select.mut[,'KRAS.mut']),'#f1eef6','#980043')
col.mat<-cbind(EGFR_col,EGFR.exon19del_col,KRAS_col,KRAS.ms12_col,KRAS.ms13_col,KRAS.ms61_col,STK11_col,KRAS_STK11_col)
colnames(col.mat)<-gsub('_col','',fixed = T,colnames(col.mat))
heatmap3(dat[complete.cases(dat[,1:7]),1:7],
dist=dist,Colv = NA,col = colorRampPalette(c('#f1eef6','red'))(1000),method = 'ward.D2',balanceColor = F,scale='none',
RowSideColors = col.mat,labCol = paste('m+',0:6),labRow = NA)
legend('topleft',lty=1,lwd=3,legend=c('co-mutation','gene-specific mutation','site-specific mutation','others'),
col=c('#980043','#dd1c77','#df65b0','#f1eef6'),cex = 0.7,bty='n',inset = c(0.1,-0.2),xpd = T)
4B-E
ecdf_mut<-function(mut,col_choice,mut_name,rest){
k<-ks.test(dat$CitG6m0[setdiff(1:nrow(dat),mut)],dat$CitG6m0[mut])
plot(ecdf(dat$CitG6m0[mut]), col=c('#df65b0','#980043')[col_choice], main=NA,
xlim=range(dat$CitG6m0,na.rm = T),xlab='CitG6m0',ylab='Cumulative Probability')
plot(ecdf(dat$CitG6m0[setdiff(1:nrow(dat),mut)]), col='#ccc6d6', add=T,main=NA)
legend(34,0.3,pch=19,col=c(c('#df65b0','#980043')[col_choice],'#ccc6d6'),legend = c(mut_name,rest),bty='n')
legend(8,0.99,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
}
par(mfrow=c(4,1),mar=c(3,3,1,1),mgp=c(2,1,0))
ecdf_mut(mut = EGFR.exon19del,col_choice = 1,mut_name = 'EGFR.exon19del',rest='Others')
ecdf_mut(mut = KRAS.ms13,col_choice = 1,mut_name = 'KRAS.ms13',rest='Others')
ecdf_mut(mut = KRAS.ms61,col_choice = 1,mut_name = 'KRAS.ms61',rest='Others')
ecdf_mut(mut = KRAS_STK11,col_choice = 2,mut_name = 'KRAS_STK11',rest='Others')
x<-dat$CitG6m0[match(c('A549','HCC44','H460'),rownames(dat))]
y<-(1/length(KRAS_STK11))*match(x,sort(dat$CitG6m0[KRAS_STK11]))
text(x+2,y,c('A549','HCC44','H460'),pos = 4,col='#980043')
4F
produce.mut<-function(x,gene){
y<-x
y[grepl(gene,x)]<-'Mut'
y[grepl('ND',x)]<-'WT'
y[grepl('NT',x)]<-NA
y[!y%in%c('Mut','WT',NA)]<-'WT'
y
}
mut$EGFR<-produce.mut(mut$Mutation,'EGFR')
dat.tumor<-dat.invivo[!grepl('adjacent',dat.invivo$Tissue.fragment,ignore.case = T),]
dat.normal<-dat.invivo[grepl('adjacent',dat.invivo$Tissue.fragment,ignore.case = T),]
dat.tumor[,3:9]<-dat.tumor[,3:9]*100
dat.normal[,3:9]<-dat.normal[,3:9]*100
ecdf_mut<-function(mut,col_choice,mut_name,rest,feature){
k<-ks.test(dat[,feature][setdiff(1:nrow(dat),mut)],dat[,feature][mut])
plot(ecdf(dat[,feature][mut]), col=c('#df65b0','#980043')[col_choice], main=NA,
xlim=range(dat[,feature]),xlab=feature,ylab='Cumulative Probability')
plot(ecdf(dat[,feature][setdiff(1:nrow(dat),mut)]), col='#ccc6d6', add=T,main=NA)
legend(34,0.3,pch=19,col=c(c('#df65b0','#980043')[col_choice],'#ccc6d6'),legend = c(mut_name,rest),bty='n')
legend(8,0.98,bty='n',legend = paste0('pv =',signif(k$p.value,2)))
}
par(oma = c(6,6,4,4) + 0.1,
mar = c(0,0,1,1) + 0.1)
layout(matrix(c(1,2,3), 1, 3, byrow = TRUE),
widths=c(5,5,2))
merge.dat<-merge(dat.tumor,mut,by.x='Pt.Code',by.y='X1',all = F)
k<-ks.test(merge.dat$`M+0`[merge.dat$EGFR=='Mut'],merge.dat$`M+0`[merge.dat$EGFR=='WT'])
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='WT']), col='#ccc6d6', main='Tumor')
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']), pch=19, cex=1.5,add=T,
xlim=c(50,100),main=NA,col="#df65b0",col.01line = "#df65b0")
lines(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']),col="#df65b0",pch=19,cex=0.6,
col.points=brewer.pal(name = 'Set1',n = 6)[match(merge.dat$Pt.Code[merge.dat$EGFR=='Mut'],
unique(merge.dat$Pt.Code[merge.dat$EGFR=='Mut']))])
legend(50,1,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
merge.dat<-merge(dat.normal,mut,by.x='Pt.Code',by.y='X1',all = F)
k<-ks.test(merge.dat$`M+0`[merge.dat$EGFR=='Mut'],merge.dat$`M+0`[merge.dat$EGFR=='WT'])
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='WT']), col='#ccc6d6', main='Adjacent Lung',yaxt='n')
plot(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']), pch=19, cex=1.5,
xlim=c(50,100),add=T,main=NA,col="#df65b0",col.01line = "#df65b0")
lines(ecdf(merge.dat$`M+0`[merge.dat$EGFR=='Mut']),col="#df65b0",pch=19,cex=0.6,
col.points=brewer.pal(name = 'Set1',n = 6)[match(merge.dat$Pt.Code[merge.dat$EGFR=='Mut'],
unique(merge.dat$Pt.Code[merge.dat$EGFR=='Mut']))])
legend(68,1,bty='n',legend = paste0('pv = ',signif(k$p.value,2)))
title(ylab = "Cumulative Probability",
outer = TRUE, line = 2.5,cex.lab=1.5)
title(xlab = "citrate m+0 ",
outer = TRUE, line = 2.5,cex.lab=1.5)
plot.new()
legend('center',pch=19,col=c('#df65b0','#ccc6d6'),legend = c('Mutant','WT'),bty='n',title = 'EGFR',xpd = T)
4H
cell_line<-sapply(lkb1$X,FUN = function(x) unlist(strsplit(x,' ',fixed = T))[1])
treatment<-sapply(lkb1$X,FUN = function(x) unlist(strsplit(x,' ',fixed = T))[2])
pv<-signif(sapply(unique(cell_line),FUN=function(cell) t.test(lkb1$M.0[cell_line==cell & treatment=='CTL'],lkb1$M.0[cell_line==cell & treatment=='LKB'])$p.value),2)
lkb1<-rbind(apply(lkb1[,-1],2,FUN=function(x) by(x,INDICES = lkb1[,1],mean)),
apply(lkb1[,-1],2,FUN=function(x) by(x,INDICES = lkb1[,1],std.error)))
lkb1<-data.frame(X=rownames(lkb1),lkb1,X.1=rep(c('mean','SEM'),each=6))
lkb1.mean<-lkb1[lkb1$X.1=='mean',2:8]
lkb1.sem<-lkb1[lkb1$X.1=='SEM',2:8]
rownames(lkb1.mean)<-rownames(lkb1.sem)<-lkb1$X[1:(nrow(lkb1)/2)]
type<-rep('',nrow(lkb1.mean))
type[grep('CTL',rownames(lkb1.mean))]<-'CTL'
type[grep('OE',rownames(lkb1.mean))]<-'LKB1 OE'
lkb1.mean$Type<-factor(type)
lkb1.mean$`CellLine`<-factor(sapply(rownames(lkb1.mean),FUN=function(x) unlist(strsplit(x,split=' ',fixed = T))[1]))
tmp<-cbind(melt(lkb1.mean),melt(lkb1.sem))[,c(1:4,6)]
colnames(tmp)[4:5]<-c('mean','sem')
dat_text <- data.frame(
x = c(1,1,1),
y = c(57,57,57),
label = pv,
CellLine=levels(tmp$CellLine)
)
tmp$variable<-gsub('M.','m+',tmp$variable,fixed = T)
ggplot(tmp,aes(x=variable,y=mean,fill=Type))+ylab('% Enrichment')+
geom_bar(stat="identity", position=position_dodge(), colour="black")+facet_wrap(~CellLine,ncol=1) +
geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem), width=.2,
position=position_dodge(.9))+theme_light()+
scale_fill_manual("Cell Line", values = c("CTL" = "black", "LKB1 OE" = "white"))+
annotate("rect", xmin = 0.75, xmax = 1.25, ymin = 50, ymax =50, alpha=1,colour = "black")+
xlab(expression('citrate from [U-'^{13}*'C] glucose labeling'))+
geom_text(data = dat_text,mapping = aes(x = x, y = y, label = label),inherit.aes = FALSE,size=2.5)
Fig 5
5C
AUC.erlotinib<-t(AUC[which(AUC$Common_chemical_name=='Erlotinib'),-(1:7)])
rownames(AUC.erlotinib)<-gsub('Calu','Calu-',rownames(AUC.erlotinib))
tmp<-merge(dat[,c('CitG6m0','CitQ6m5')],
S4.data[,c('EGFR-pY1173','E-cadherin','Beta-catenin')],
by.x=0,by.y = 0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
tmp<-merge(tmp,AUC.erlotinib,
by.x=0,by.y = 0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
colnames(tmp)[6]<-'Erlotinib AUC'
tmp.5C<-tmp
tmp.col<-rep('black',nrow(tmp))
tmp.col[rownames(tmp)%in%c('H1869','H1650','HCC2935','HCC1897','HCC4019')]<-'purple'
tmp.col[rownames(tmp)%in%c('H1299','H1355','DFCI032','H1975','H23')]<-'green'
tmp$col<-factor(tmp.col)
colnames(tmp)<-gsub(' ','_',colnames(tmp),fixed = T)
colnames(tmp)<-gsub('-','_',colnames(tmp),fixed = T)
g<-ggpairs(tmp,
upper=list(continuous=cor_fun),
lower=list(mapping = aes(color=col,alpha=0.8),size=0.05),
columns = 1:6)+theme_bw()+theme(panel.grid = element_blank())
colors<-c('darkgray','#4d9221','#c51b7d')
for(i in 2:6){
for(j in 1:(i-1)){
g[i,j]<-g[i,j] + scale_color_manual(values=colors)
}
}
g
5E-figure
col.grad<-function(cols=c('blue','gray','red'),x){
y=scales::alpha(colorRampPalette(cols)(length(x)),alpha = 1)[rank(x)]
y[is.na(x)]<-'gray'
return(y)
}
color.bar <- function(lut, min, max=-min, ticks=seq(min, max, len=2), title='') {
scale = (length(lut)-1)/(max-min)
plot(c(min,max), c(0,1), type='n', bty='n', xaxt='n', xlab='', yaxt='n', ylab='', main=title,mar=c(0,3,2,3),ylim=c(0,4))
axis(3, ticks, las=1,labels = c('low','high'))
for (i in 1:(length(lut)-1)) {
x = (i-1)/scale + min
rect(x,3,x+1/scale,4, col=lut[i], border=NA)
}
}
tmp<-tmp.5C
colnames(tmp)[6]<-'AUC'
tmp$Name<-rownames(tmp)
tmp.h<-tmp[rownames(tmp)%in%c('H1869','H1650','HCC2935','HCC1897','HCC4019'),]
tmp.l<-tmp[rownames(tmp)%in%c('H1299','H1355','DFCI032','H1975','H23'),]
mut.shape=ifelse(mutations[match(rownames(tmp),rownames(mutations)),'EGFR']=="",19,2)
mut.shape[!is.na(str_extract(pattern = "p.[A-Z]7[0-9]+_[A-Z]7[0-9]+del",
string = mutations[match(rownames(tmp),rownames(mutations)),'EGFR']))]<-17
AUC.size<-rep(3,nrow(tmp))
AUC.size[is.na(AUC.erlotinib[match(rownames(tmp),rownames(AUC.erlotinib))])]<-1.5
ggplot(tmp, aes(x= CitG6m0, y = CitQ6m5, label=Name)) +
geom_point(
colour=col.grad(cols=rev(brewer.pal(7,'PuOr')),
x=tmp$AUC),
size=AUC.size,
shape=mut.shape
) +
geom_text_repel(data = tmp.h,
colour='#c51b7d',
nudge_x = -tmp.h$CitG6m0,
segment.size = 0.2,
segment.color = "grey50") +
geom_text_repel(data = tmp.l,
colour="#4d9221",
nudge_x = 70-tmp.l$CitG6m0,
segment.size = 0.2,
segment.color = "grey50") +
theme_classic(base_size = 15)
5E-legend
color.bar(colorRampPalette(rev(brewer.pal(7,'PuOr')))(100), -1)
legend('bottomleft',pch=c(19,17,2),legend = c('WT','exon19.del','Other'),title = 'EGFR status',bty = 'n')
legend('bottomright',pch=19,pt.cex = c(1,0.5),legend = c("available","missing"),col=c('black','grey'),title = 'Erlotinib AUC',bty = 'n')
5D
ge.emt<-ge[intersect(rownames(dat)[complete.cases(dat$CitG6m0)],rownames(ge)),intersect(colnames(ge),emt.signature)]
ge.emt[]<-apply(ge.emt,2,scale)
col.mat<-apply(dat[rownames(ge.emt),c('CitG6m0','CitQ6m5')],2,
FUN=function(x) col.grad(cols=rev(brewer.pal(7,'PiYG')),x=x))
heatmap3(ge.emt,
RowSideColors = col.mat,
dist=dist,method = 'ward.D2',labCol = NA,labRow = NA)
5I-preparation and model fitting
tmp<-tmp.5C
colnames(tmp)[6]<-'AUC'
emt.pca<-prcomp(ge.emt)
mc<-Mclust(emt.pca$x[,1],G=2)$classification
tmp$`EMT-Signature`<-mc[match(rownames(tmp),rownames(ge.emt))]
tmp<-tmp[complete.cases(tmp),]
fit<-lm(AUC~.,tmp)
step<-stepAIC(fit, direction="both")
## Start: AIC=446.85
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `Beta-catenin` +
## `EMT-Signature`
##
## Df Sum of Sq RSS AIC
## - `Beta-catenin` 1 1212 215542 445.14
## - CitQ6m5 1 1911 216240 445.31
## - `EMT-Signature` 1 2807 217137 445.53
## <none> 214329 446.85
## - `E-cadherin` 1 10665 224995 447.37
## - CitG6m0 1 16888 231218 448.79
## - `EGFR-pY1173` 1 50876 265205 455.92
##
## Step: AIC=445.14
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `EMT-Signature`
##
## Df Sum of Sq RSS AIC
## - CitQ6m5 1 1639 217181 443.54
## - `EMT-Signature` 1 2293 217835 443.69
## <none> 215542 445.14
## - `E-cadherin` 1 9883 225425 445.47
## - CitG6m0 1 15682 231224 446.79
## + `Beta-catenin` 1 1212 214329 446.85
## - `EGFR-pY1173` 1 50231 265773 454.04
##
## Step: AIC=443.54
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin` + `EMT-Signature`
##
## Df Sum of Sq RSS AIC
## - `EMT-Signature` 1 2890 220071 442.22
## <none> 217181 443.54
## - `E-cadherin` 1 11954 229135 444.32
## + CitQ6m5 1 1639 215542 445.14
## + `Beta-catenin` 1 941 216240 445.31
## - CitG6m0 1 21736 238917 446.50
## - `EGFR-pY1173` 1 48735 265915 452.06
##
## Step: AIC=442.22
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin`
##
## Df Sum of Sq RSS AIC
## <none> 220071 442.22
## - `E-cadherin` 1 11723 231794 442.92
## + `EMT-Signature` 1 2890 217181 443.54
## + CitQ6m5 1 2236 217835 443.69
## + `Beta-catenin` 1 424 219647 444.12
## - CitG6m0 1 19444 239515 444.63
## - `EGFR-pY1173` 1 50407 270478 450.95
5I-stepwise feature selection result
step$anova
## Stepwise Model Path
## Analysis of Deviance Table
##
## Initial Model:
## AUC ~ CitG6m0 + CitQ6m5 + `EGFR-pY1173` + `E-cadherin` + `Beta-catenin` +
## `EMT-Signature`
##
## Final Model:
## AUC ~ CitG6m0 + `EGFR-pY1173` + `E-cadherin`
##
##
## Step Df Deviance Resid. Df Resid. Dev AIC
## 1 45 214329.3 446.8493
## 2 - `Beta-catenin` 1 1212.298 46 215541.6 445.1426
## 3 - CitQ6m5 1 1639.044 47 217180.7 443.5365
## 4 - `EMT-Signature` 1 2890.472 48 220071.1 442.2240
5I-table
fit1<-lm(AUC~`CitG6m0`+`EGFR-pY1173`+`E-cadherin`,data = tmp)
fit2<-lm(AUC~`CitG6m0`+`EMT-Signature`+`EGFR-pY1173`+`E-cadherin`,data=tmp)
tab_model(fit1,fit2,
dv.labels = c('Model 1 (AIC selected)','Model 2 (added EMT-Signature)'))
|
Â
|
Model 1 (AIC selected)
|
Model 2 (added EMT-Signature)
|
|
Predictors
|
Estimates
|
CI
|
p
|
Estimates
|
CI
|
p
|
|
(Intercept)
|
553.43
|
503.88 – 602.98
|
<0.001
|
525.98
|
441.71 – 610.25
|
<0.001
|
|
Cit G 6 m 0
|
-1.89
|
-3.69 – -0.09
|
0.045
|
-2.04
|
-3.89 – -0.20
|
0.035
|
EGFR-pY1173
|
-28.06
|
-44.64 – -11.47
|
0.002
|
-27.64
|
-44.32 – -10.96
|
0.002
|
E-cadherin
|
-11.22
|
-24.97 – 2.53
|
0.116
|
-18.26
|
-40.52 – 3.99
|
0.114
|
EMT-Signature
|
|
|
|
28.71
|
-42.43 – 99.84
|
0.433
|
|
Observations
|
52
|
52
|
|
R2 / adjusted R2
|
0.457 / 0.423
|
0.464 / 0.418
|
5J
yy<-data.frame(x=fit1$fitted.values,y=tmp$AUC)
ggplot(yy,aes(x=x,y=y))+geom_point(colour='navy',alpha=0.4,size=3)+
geom_smooth(method='lm',se=F,colour=scales::alpha('black',0.3))+
xlab('fitted values')+ylab('measured values')+ggtitle('Erlotinib AUC')+theme_classic()
Fig 6
6A
c<-rcorr(as.matrix(dat[,c('SerG6m3','GlyG6m2')]))
ggplot(dat[,c('SerG6m3','GlyG6m2')],mapping = aes(x=SerG6m3,y=GlyG6m2))+
geom_point(colour='#6e016b',alpha=0.5,size=2.5)+
theme_bw()+
ggtitle(paste0("Serine Glycine interconversion\nr = ",round(c$r[1,2],2),', pv = ',signif(c$P[1,2],2)))
6B
drug<-S4.data[,grep(")$",colnames(S4.data))]
rownames(drug)<-S4.data$Cell.Line
drug<-as.matrix(drug[match(rownames(dat),rownames(drug)),]);class(drug)<-'numeric'
pv<-apply(drug,2,FUN=function(x) cor.test(dat$SerG6m3,x,method = 'pearson')$p.value)
r<-apply(drug,2,FUN=function(x) cor.test(dat$SerG6m3,x,method = 'pearson')$estimate)
n<-apply(drug,2,FUN=function(x) length(which(!is.na(x) & !is.na(dat$SerG6m3))))
res<-data.frame(pv,r,n)
rownames(res)<-sapply(rownames(res),FUN=function(x) unlist(strsplit(x,split = '.',fixed = T))[1])
res<-res[order(res$pv,decreasing = F),]
res$pv<-(-log10(res$pv))
res$drug<-factor(rownames(res),levels = rev(rownames(res)))
res$significance<-rep('insig',nrow(res))
res$significance[1:countSignificant(Bum(pv),0.05)]<-'sig'
res$significance<-factor(res$significance,levels = c('insig','sig'))
ggplot(res,aes(y=pv,fill=significance,x=drug))+geom_bar(stat="identity")+coord_flip() +scale_fill_brewer(palette="PuRd")+
geom_hline(yintercept = (-log10(0.05)), linetype="dotted")+theme_classic()+ylab('-log10(p-value)')+xlab("")+ggtitle('Correlation with SerG6m3')
6D
tmp<-data.frame(Pemetrexed=drug[,"Pemetrexed.(uM)"],SerG6m3=dat$SerG6m3)
rownames(tmp)<-rownames(dat)
tmp<-tmp[complete.cases(tmp),]
mc<-Mclust(tmp$Pemetrexed,G=2)
tmp$sensitivity<-c('sensitive','resistant')[mc$classification]
tmp$sensitivity<-factor(tmp$sensitivity,levels=c('sensitive','resistant'))
tmp$cutoff<-c('low','high')[as.numeric(tmp$SerG6m3>16)+1]
h.cells<-c('PC-9','H2170','H596','H460','H1299')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195')
tmp$name<-rownames(tmp)
tmp.h<-tmp[rownames(tmp)%in%h.cells,]
tmp.l<-tmp[rownames(tmp)%in%l.cells,]
h.col<-brewer.pal(n = 3,'Set1')[1]
l.col<-brewer.pal(n = 3,'Set1')[2]
ggplot(tmp,aes(x=SerG6m3,y=Pemetrexed,colour=sensitivity,shape=cutoff, label=name))+geom_point()+
geom_text_repel(data = tmp.h,
colour=h.col,
nudge_y = 1-tmp.h$Pemetrexed,
segment.size = 0.2,
segment.color = "grey50") +
geom_text_repel(data = tmp.l,
colour=l.col,
nudge_x = 70-tmp.l$SerG6m3,
segment.size = 0.2,
segment.color = "grey50") +
theme_bw()
6E
mat.mean<-function(x) apply(x,2,FUN=function(y) mean(y,na.rm = T))
pem.mean<-do.call(rbind,by(data = pem[,-1],INDICES = pem[,1],FUN = mat.mean))
pem.auc<-apply(pem.mean,2,FUN=function(x) trapz(log10(as.numeric(rownames(pem.mean))),x))
pem.plot<-data.frame(cell=names(pem.auc),auc=pem.auc)
h.cells<-c('PC-9','H2170','H596','H460','H1299')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195')
pem.plot$SerG6m3<-rep(NA,nrow(pem.plot))
pem.plot$SerG6m3[as.character(pem.plot$cell)%in%h.cells]<-'high'
pem.plot$SerG6m3[as.character(pem.plot$cell)%in%l.cells]<-'low'
pem.plot$SerG6m3<-factor(pem.plot$SerG6m3,levels = c('high','low'))
pem.plot$cell<-factor(pem.plot$cell,levels = c(h.cells,l.cells))
pem.plot$mean<-unlist(by(data = pem.plot$auc,INDICES = pem.plot$SerG6m3,mean))[match(pem.plot$SerG6m3,c('high','low'))]
pem.plot$sd<-unlist(by(data = pem.plot$auc,INDICES = pem.plot$SerG6m3,sd))[match(pem.plot$SerG6m3,c('high','low'))]
set.seed(621)
ggplot(pem.plot,aes(x=SerG6m3,y=auc,colour=SerG6m3))+ geom_jitter(width = 0.2)+scale_color_brewer(palette="Set1")+
ggtitle(paste0('pv = ',signif(t.test(auc~SerG6m3,pem.plot)$p.value,1)))+ylab('AUC')+
geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd, group = SerG6m3),width = 0.2)+theme_bw()
6F-H
g.list<-list()
for(sheet in sheets(wb.misty)){
tmp<-readWorkbook(wb.misty,sheet,rowNames = T)
tmp<-data.frame(days=rownames(tmp),Cell=sheet,melt(tmp))
tmp$variable<-sapply(as.character(tmp$variable),FUN=function(x) substr(x,1,nchar(x)-2))
tmp$variable<-factor(tmp$variable,levels = c('Control','Pem.500', 'Pem.750','Pem.1000'))
tmp<-tmp[complete.cases(tmp),]
tmp$days<-as.character(tmp$days)
colnames(tmp)[colnames(tmp)=='variable']<-'Treatment'
g.list[[sheet]]<-ggpar(ggline(tmp[tmp$Cell==sheet,], x = "days", y = "value", add = "mean_se",color = 'Treatment',
palette = "Dark2")+ylab('tumor volume')+
stat_compare_means(aes(group = Treatment), label = "p.signif", method = 'aov',
label.y = rep(max(na.omit(unlist(by(data = tmp$value,INDICES = list(tmp$Treatment,tmp$days),mean))))*1.2,
length(unique(tmp$days[tmp$Cell==sheet]))))+
ggtitle(sheet),legend='none')
}
marrangeGrob(g.list, nrow=1, ncol=3,top=NULL)
6F-H-legend
legend <- cowplot::get_legend(ggpar(ggline(tmp[tmp$Cell==sheet,], x = "days", y = "value", add = "mean_se",color = 'Treatment',
palette = "Dark2")+ylab('tumor volume')+
stat_compare_means(aes(group = Treatment), label = "p.signif", method = 'aov',
label.y = rep(max(na.omit(unlist(by(data = tmp$value,INDICES = list(tmp$Treatment,tmp$days),mean))))*1.2,
length(unique(tmp$days[tmp$Cell==sheet]))))+
ggtitle(sheet),legend='right'))
grid.newpage();grid.draw(legend)
Fig S1
S1A
doubling.time=(48/log2(dat$`Day3/Day1`)+96/log2(dat$`Day5/Day1`))/2
nutrients<-dat[,c('Lac','Glc','Lac/Glc','Glu','Gln','Glu/Gln')]
nutrients.norm<-nutrients
nutrients.norm[]<-apply(nutrients,2,FUN=function(x) x/((2^(7/doubling.time)+1)/2))
nutrients.norm$`Glu/Gln`<-nutrients.norm$Glu/nutrients.norm$Gln
nutrients.norm$`Lac/Glc`<-nutrients.norm$Lac/nutrients.norm$Glc
tmp<-data.frame(melt(nutrients[,c(1,2,4,5)]),corrected=melt(nutrients.norm[,c(1,2,4,5)])[,2])
colnames(tmp)[1:2]<-c('feature','uncorrected')
tmp$cor<-rep(unlist(lapply(names(nutrients)[c(1,2,4,5)],FUN=function(x) round(cor.test(nutrients[,x],nutrients.norm[,x])$estimate,3))),each=nrow(nutrients))
tmp$feature_name=factor(paste(tmp$feature,tmp$cor,sep="\nr = "),ordered = F)
tmp$feature_name=factor(tmp$feature_name,levels(tmp$feature_name)[c(1,4,2,3)])
ggplot(tmp,aes(x=uncorrected,y=corrected))+geom_point(alpha=0.3)+theme_classic()+facet_wrap(~feature_name,scale='free',ncol = 1)+
xlab('normalized to final protein content')+ylab('normalized to estimated average protein content in 7h growth')
S1B
replicate.sd<-list()
cell.sd<-list()
new.names<-c('CitG6','CitG24','CitQ6','CitQ24','SerG6','SerG24',
'GlyG6','GlyG24','FumG6','FumG24','FumQ6','FumQ24',
'MalG6','MalG24','MalQ6','MalQ24','LacG6','LacG24','LacQ6','LacQ24')
result.list<-list()
sd.mat.list<-list()
for(i in 2:(length(new.names)+1)){
sheet<-sheets(wb)[i]
tmp<-read.xlsx(wb, sheet = sheet,rowNames = T,colNames = T)
sd.mat<-as.matrix(tmp[-1,grep('Std',colnames(tmp),ignore.case = T):(grep('Rep',colnames(tmp),ignore.case = T)-1)])
rownames(sd.mat)<-tmp[-1,1]
sd.mat<-sd.mat[complete.cases(sd.mat),]
class(sd.mat)<-'numeric'
btw.cell.sd<-apply(as.matrix(tmp[-1,grep('ave',colnames(tmp),ignore.case = T):(grep('Std',colnames(tmp),ignore.case = T)-1)]),2,FUN=function(x) sd(x,na.rm = T))
result.list[[length(result.list)+1]]<-data.frame(feature=rep(sheet,length(btw.cell.sd)),btw.cell.sd,btw.rep.sd.mean=apply(sd.mat,2,mean),btw.rep.sd.sd=apply(sd.mat,2,sd))
sd.mat.list[[length(sd.mat.list)+1]]<-sd.mat
names(sd.mat.list)[length(sd.mat.list)]<-new.names[i-1]
}
tmp<-do.call(rbind,result.list)
tmp$metabolite<-substr(tmp$feature,1,3)
tmp$`labeling time`<-factor(paste(substr(tmp$feature,5,8),'hr'),levels=c('6 hr','24 hr'))
tmp$tracer<-c('glucose','glutamine')[match(substr(tmp$feature,4,4),c('G','Q'))]
ggplot(tmp,aes(x=btw.cell.sd,y=btw.rep.sd.mean,colour=metabolite,shape=tracer,alpha=`labeling time`))+
geom_abline(aes(intercept=0,slope=1),colour='lightgray')+
scale_alpha_discrete(range=c(0.3,0.8))+
geom_point()+coord_flip()+
geom_pointrange(aes(ymin=btw.rep.sd.mean, ymax=btw.rep.sd.mean+btw.rep.sd.sd))+
theme_classic()+ylim(0,22)+xlim(0,22)
S1C was made by Pei-Hsuan Chen
S1D
diversity.plot<-function(sheet,cell_lines,tracer){
tmp<-read.xlsx(wb, sheet = sheet,rowNames = T,colNames = T,startRow = 2)
tmp$Cell.line<-toupper(tmp$Cell.line)
tmp<-tmp[tmp$Cell.line%in%cell_lines,]
rownames(tmp)<-tmp$Cell.line;tmp<-tmp[,-c(1:ifelse(sheet=='CitG6',16,15))]
tmp.mean<-apply(t(tmp),2,FUN=function(x){
apply(matrix(byrow = F,data = x,nrow = 7,ncol = length(x/7)),1,FUN = function(y) mean(na.omit(y)))
})
tmp.sem<-apply(t(tmp),2,FUN=function(x){
apply(matrix(byrow = F,data = x,nrow = 7,ncol = length(x/7)),1,FUN = function(y) std.error(na.omit(y)))
})
rownames(tmp.sem)<-rownames(tmp.mean)<-paste('m+',0:6,sep = '')
tmp.m<-setNames(cbind(melt(tmp.mean),melt(tmp.sem)[,3]),c('isotopologue','cell line','mean','sem'))
if(tracer=='glucose'){
title<-expression('citrate labeling by [U-'^{13}*'C] glucose labeling')
}else{
title<-expression('citrate labeling by [U-'^{13}*'C] glutamine labeling')
}
ggplot(tmp.m,aes(x=isotopologue,y=mean))+
geom_bar(stat="identity",position=position_dodge())+
geom_errorbar(aes(ymin=mean-sem,ymax=mean+sem),width=.2,position=position_dodge(.9))+
ylab('% enrichment')+
ggtitle(title)+
theme_bw()+
theme(axis.text.x = element_text(angle = 45, hjust = 1))+
facet_grid(.~`cell line`)
}
diversity.plot(sheet='CitG6',cell_lines = c('HCC44','H1869','H1395'),tracer='glucose')
S1E
diversity.plot(sheet='CitQ6',cell_lines = c('H1975','H441','H1395'),tracer='glutamine')
Fig S5
S5C
h<-hclust(dist(ge.emt),method = 'ward.D2')
hc<-cutree(h,2)
# the line commented below was used to check which cluster was the epithelial cluster
# boxplot(ge.emt[,'CDH1']~hc)
E<-rownames(ge.emt)[hc==2]
M<-rownames(ge.emt)[hc==1]
metab.p<-unlist(lapply(c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln'),
FUN=function(x) t.test(na.omit(dat[rownames(dat)%in%E,x]),na.omit(dat[rownames(dat)%in%M,x]))$p.value))
names(metab.p)<-c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')
tmp<-data.frame(dat[match(c(E,M),rownames(dat)),c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')],
emt=rep(c('epithelial','mesenchymal'),c(length(E),length(M))),check.names = F)
dat_text <- data.frame(
x = rep(1.85,4),
y = 0.9*apply(tmp[,1:4],2,FUN=function(x) max(x,na.rm = T)),
label = paste("pv =",signif(metab.p,2)),
variable=c("CitG6m0","CitQ6m5",'Lac/Glc','Glu/Gln')
)
tmp<-melt(tmp)
ggplot(tmp,aes(x=emt,y=value))+geom_boxplot()+
geom_text(data = dat_text,mapping = aes(x = x, y = y, label = label),inherit.aes = FALSE,size=4)+
facet_wrap(~variable,scale='free_y')+xlab("")+theme_bw()
S5D
tmp<-merge(CDH1.data,dat[,'CitG6m0',drop=F],by.x=0,by.y=0,all = T)
rownames(tmp)<-tmp[,1];tmp<-tmp[,-1]
tmp<-tmp[complete.cases(tmp),]
ggpairs(tmp[complete.cases(tmp),],
upper=list(continuous=cor_fun),
lower=list(mapping = aes(alpha=0.8),size=0.05))+
theme_bw()+theme(panel.grid = element_blank())
S6A
common<-intersect(rownames(dat)[!is.na(dat$SerG6m3)],rownames(ge))
dat.Ser<-dat[common,]
ge.Ser<-ge[common,c('PHGDH','PSAT1','PSPH','CBS','SHMT1','SHMT2')]
ge.Ser<-ge.Ser[order(dat.Ser$SerG6m3,decreasing = T),]
ge.Ser[]<-apply(ge.Ser,2,rank)
annot_df<-dat.Ser[order(dat.Ser$SerG6m3,decreasing = T),'SerG6m3',drop=F]
annot_df[,1]<-rank(annot_df[,1])
mf_range<-range(annot_df[,1])
col<-circlize::colorRamp2(seq(from=mf_range[1],to=mf_range[2],length.out = 3),
c(l.col, "#EEEEEE", h.col))
col<-list(col)
names(col)<-'SerG6m3'
# Create the heatmap annotation
ha <- HeatmapAnnotation(annot_df, col = col, show_legend = F)
colnames(ge.Ser)<-paste(colnames(ge.Ser),
"\n",'rho=',round(apply(ge.Ser,2,FUN=function(x) cor.test(x,sort(dat.Ser$SerG6m3,decreasing = T),method = 'spearman')$estimate),1),
', p=',signif(apply(ge.Ser,2,FUN=function(x) cor.test(x,sort(dat.Ser$SerG6m3,decreasing = T),method = 'spearman')$p.value),1),sep = '')
ht<-Heatmap(t(ge.Ser),
show_row_names = T,show_column_names = F,
cluster_rows = F,cluster_columns = F,
show_heatmap_legend = F,
top_annotation = ha)
lgd = list(Legend(at = c("low","","medium","","high"), title = "gene expr", type = "grid",
legend_gp = gpar(fill = c("#0000FF","#7777F6","#EEEEEE","#F67777","#FF0000"))),
Legend(at = c("low",'medium',"high"), title = "SerG6m3", type = "grid",
legend_gp = gpar(fill = c(l.col, "#EEEEEE", h.col))))
draw(ht, column_title = 'Serine metabolism genes',annotation_legend_list = lgd)
an<-'SerG6m3'
decorate_annotation(an, {
# annotation names on the left
grid.text(an, unit(1, "npc") + unit(2, "mm"), 0.5, default.units = "npc", just = "left")
})
S6B
tmp.cor<-rcorr(as.matrix(dat[,c('SerG6m3','Day3/Day1')]))
ggplot(data = dat[,c('SerG6m3','Day3/Day1')],aes(x=SerG6m3,y=`Day3/Day1`))+geom_point(color=scales::alpha('#8856a7',0.7))+
ggtitle(paste0('r = ',round(tmp.cor$r[1,2],2),', pv = ',signif(tmp.cor$P[1,2],2)))+theme_bw()
S6C-J
h.cells<-c('PC-9','H2170','H596','H460','H1299','H1155')
l.cells<-c('HCC44','H2882','H650','H1395','HCC1195','H2250')
hl.plot<-function(tmp,feature,hide=T){
tmp.plt<-data.frame(cell=names(tmp),value=tmp)
tmp.plt$SerG6m3<-rep(NA,nrow(tmp.plt))
tmp.plt$SerG6m3[tmp.plt$cell%in%h.cells]<-'high'
tmp.plt$SerG6m3[tmp.plt$cell%in%l.cells]<-'low'
tmp.plt$SerG6m3<-factor(tmp.plt$SerG6m3,levels = c('high','low'))
tmp.plt$mean<-unlist(by(data = tmp.plt$value,INDICES = tmp.plt$SerG6m3,mean))[match(tmp.plt$SerG6m3,c('high','low'))]
tmp.plt$sd<-unlist(by(data = tmp.plt$value,INDICES = tmp.plt$SerG6m3,sd))[match(tmp.plt$SerG6m3,c('high','low'))]
tmp.plt<-tmp.plt[complete.cases(tmp.plt),]
pv<-t.test(value~SerG6m3,tmp.plt)$p.value
g<-ggplot(tmp.plt,aes(x=SerG6m3,y=value,colour=SerG6m3))+ geom_jitter(width = 0.2)+scale_color_brewer(palette="Set1")+
ggtitle(paste('pv =',signif(pv,2)))+ylab(feature)+
geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd, group = SerG6m3),width = 0.2)+theme_bw()
if(hide) g+theme(legend.position = "none")
else g
}
g.list<-list()
g.list[[1]]<-hl.plot(apply(brdu,2,mean),feature='Brdu Staining')
for(i in 1:4){
g.list[[i+1]]<-hl.plot(apply(as.matrix(nova[[i]]),2,mean),feature=names(nova)[i])
}
g.list[[6]]<-hl.plot(apply(ser.pool,2,mean),feature='Serine Pool Size')
g.list[[7]]<-hl.plot(apply(as.matrix(pt.all$Ser),2,mean),feature='Intracellular Serine')
g.list[[8]]<-hl.plot(apply(as.matrix(med.all$Ser),2,mean),feature='Extracellular Serine')
marrangeGrob(g.list, layout_matrix = matrix(1:9, byrow = T, nrow = 3),top=NULL)
S6C-J-legend
legend <- cowplot::get_legend(hl.plot(apply(brdu,2,mean),feature='Brdu Staining',hide = F))
grid.newpage();grid.draw(legend)
